function sys = phase_peak(domain, fn, bw, phase, type)

if strcmp(domain{1}, "s")
    
elseif strcmp(domain{1}, "z")
    ts = domain{2};
    fs = 1/ts;
end

phase_max = phase / 180 * pi;
phase_max_omega = fn * 2 * pi;
peak_notch_bw = bw * 2 *pi;

phase_max_omega_z = phase_max_omega / fs;
phase_max_omega = 2 * fs * tan(phase_max_omega_z / 2);
peak_notch_bw_z = peak_notch_bw / fs;
peak_notch_bw = 2 * fs * tan(peak_notch_bw_z / 2);

a1(1) = 1;
a1(2) = peak_notch_bw;
a1(3) = -phase_max_omega^2;
tmp1 = roots(a1);
amp_min_omega = max(tmp1);
amp_max_omega = amp_min_omega + peak_notch_bw;

x = sqrt((((tan((pi - phase_max) / 2))^2) * (peak_notch_bw^2)) / (4 * phase_max_omega^2 + ((tan((pi - phase_max) / 2))^2) * 2 * (peak_notch_bw^2 + 2 * phase_max_omega^2)));



a2(1) = 1;
a2(2) = -(1 - 2 * x^2) * (amp_min_omega^2 + amp_max_omega^2);
a2(3) = amp_min_omega^2 * amp_max_omega^2;
tmp2 = roots(a2);
wz = sqrt(min(tmp2));

wp = amp_min_omega * amp_max_omega / wz;

xz = x;
xp = x;

Cs = tf(wp^2 / wz^2 * [1 2*xz*wz wz^2],[1 2*xp*wp wp^2]);
if type == 1
    Cs = Cs / wp^2 * wz^2;
end

if phase <= 0
    Cs = 1/Cs;
end

Cz = c2d(Cs, 1/fs, 't');
[num_z, den_z] = tfdata(Cz, 'v');


if strcmp(domain{1}, "s")
    sys = Cs;
elseif strcmp(domain{1}, "z")
    sys = Cz;
end

end